Signatures of Cholera Outbreak during the Yemeni Civil War, 2016–2019

The Global Task Force on Cholera Control (GTFCC) created a strategy for early outbreak detection, hotspot identification, and resource mobilization coordination in response to the Yemeni cholera epidemic. This strategy requires a systematic approach for defining and classifying outbreak signatures, or the profile of an epidemic curve and its features. We used publicly available data to quantify outbreak features of the ongoing cholera epidemic in Yemen and clustered governorates using an adaptive time series methodology. We characterized outbreak signatures and identified clusters using a weekly time series of cholera rates in 20 Yemeni governorates and nationally from 4 September 2016 through 29 December 2019 as reported by the World Health Organization (WHO). We quantified critical points and periods using Kolmogorov–Zurbenko adaptive filter methodology. We assigned governorates into six clusters sharing similar outbreak signatures, according to similarities in critical points, critical periods, and the magnitude of peak rates. We identified four national outbreak waves beginning on 12 September 2016, 6 March 2017, 28 May 2018, and 28 January 2019. Among six identified clusters, we classified a core regional hotspot in Sana’a, Sana’a City, and Al-Hudaydah—the expected origin of the national outbreak. The five additional clusters differed in Wave 2 and Wave 3 peak frequency, timing, magnitude, and geographic location. As of 29 December 2019, no governorates had returned to pre-Wave 1 levels. The detected similarity in outbreak signatures suggests potentially shared environmental and human-made drivers of infection; the heterogeneity in outbreak signatures implies the potential traveling waves outwards from the core regional hotspot that could be governed by factors that deserve further investigation.


Introduction
The ongoing cholera epidemic in Yemen is the largest in recorded history, with over 1 million suspected cholera cases since 2016 [1,2]. The inability of the Yemeni government, Houthi forces, and other involved actors and international agencies to mitigate this epidemic demonstrates the challenges that many countries face when trying to manage infectious disease outbreaks amidst other humanitarian emergencies [1][2][3][4][5]. The ongoing outbreak coincides with the Yemeni Civil War (2015-ongoing), which has crippled existing health infrastructure and depleted medical resource stockpiles [5,6]. Even before these crises, Yemen was one of the poorest countries in the Arabian Peninsula, with widespread malnutrition and rampant poverty [5].

Time Series Compilation
To conduct our analyses, we compiled all abstracted records in a uniform format. We aggregated data by WHO-defined epidemiological weeks as performed for other WHOmonitored infections [32]. Daily case reports were available for 29 weeks (17% of study period) from Week  We extracted daily case data from irregularly reported WHO EMRO daily epidemiological bulletins and aggregated cases by WHO-defined weeks. Bulletins provided updates on the number of cumulative confirmed infections occurring during multi-day reporting periods based on time of hospital admission [26,27]. We estimated average daily cases by subtracting total cases from consecutive reports and dividing by the number of days in the reporting period. We prorated weekly case counts for Week 12 of 2017 as records were incomplete (5 of 7 days had missing records). We interpolated records for 4 consecutive weeks (14% out of 29) using linear approximation based on information from adjacent weeks (before and after the missing week) capturing the overall seasonal behavior.
We extracted weekly records from epidemiological bulletins available from Week 17 of 2017 (24 April 2017) through Week 26 of 2018 (1 July 2018). Weekly data were reported using a governorate-level spreadsheet containing cumulative confirmed cholera cases for the WHO-defined "Second Wave" of the epidemic. In these bulletins, weekly cases were provided for the week of the report, including both provisional and confirmed cases, and for each of the prior three weeks, which included only confirmed cases. We reviewed each weekly bulletin in chronological order and extracted the 3 week lagged confirmed case estimate that presumably contained most accurate information. For 5 incomplete bulletins (8% out of 62), we estimated weekly cases using 2 week and 1 week lagged case estimates. We interpolated cases for Week 8 of 2018 (19)(20)(21)(22)(23)(24)(25) using adjacent weeks' values.
We used monthly records from the WHO EMRO monthly situation reports available from July 2018 through December 2019. These bulletins provided cumulative confirmed cases from Week 17 of 2017 to the final day of the reporting month. We converted monthly cumulative counts to weekly case counts in two steps. First, we subtracted consecutive monthly bulletins to estimate newly reported cases in that month. Next, we estimated average cases for each day within a given month by dividing these new cases by the number of days within the month, and aggregated cases according to WHO-defined epidemiological weeks. We interpolated missing records for 4 consecutive weeks (5.6% out of 72) using linear approximation based on information from adjacent weeks.
We estimated national cases by summing all governorate-level cases for each week of the 172 week study. We conducted all analyses according to study week (ranging 1-172) and reported all results using a combination of WHO-defined epidemiological weeks (ranging 1-52 per year) and Gregorian calendar dates (formatted DD-Month-YYYY).

Rate Calculations
To compare outbreak signatures across governorates, we calculated weekly rates of infection using fatality-adjusted, pro-rated weekly population estimates. We extracted data on conflict fatalities using the Armed Conflict Location and Event Data (ACLED) project, which reported on both violent and non-violent political events [33,34]. We summed all fatalities from all political events on each day per governorate and then aggregated daily totals by WHO-defined epidemiological week. We used fatality counts to correct the population estimates. We also applied a low-to-moderate population growth rate (≈ 0.024) reported by the Yemeni Central Bureau of Statistics to approximate population changes during the civil war [35] and prorated this growth rate for the annual cycle in weeks (e.g., 0.024/52.25 = 0.00045933).
We estimated weekly governorate-level population in several steps. First, we calculated the population for Week [30,35,36]. Next, for Week 2 of 2017 through Week 52 of 2019, we estimated weekly governorate-level population iteratively forwards by subtracting conflict fatalities from the prior week and multiplying the difference by the prorated annual growth rate (Equation (1)). Finally, to correct for fatalities occurring between Week 36 of 2016 and Week 1 of 2017, we estimated weekly population by adding conflict fatalities and dividing the sum by the weekly-adjusted population growth rate (Equation (2)). P s,t = (P s,t−1 − F s,t−1 ) * (1 + 0.024 52.25 ) if t > Week 1 of 2017 (1) P s,t = (P s,t+1 + F s,t+1 ) * (1 − 0.024 52.25 ) if t < Week 1 of 2017 (2) where P s,t -the population for s-governorate in t-week; and F s,t -the all-cause war conflict fatalities for s-governorate in t-week. We calculated weekly rates of confirmed cholera infections (cases per 100,000 persons or cph) by dividing weekly cases by population estimates and multiplying by 100,000. We provide the curated weekly time series of case and rate estimates with and without interpolated values (Supplementary Excel Table S1).

Defining Outbreak Signatures and Features
We defined outbreak signatures using the timing of critical points, the duration of critical periods, and the magnitude of smoothed peak rates (example shown in Figure 1; Supplementary Excel Table S2). Critical points included the onset (O) of an outbreak when weekly rates began increasing, an outbreak peak (P) when weekly rates reached local maxima, and the return-from-peak or resolution (R) of an outbreak when weekly rates returned to steady-state levels, if this occurred.  (trivial derivative) where acceleration periods are shown by the orange gradient, deceleration periods are shown by the blue gradient, and steady-state periods or critical change points near peaks are shown in grey. All panels share the common horizontal axis of time (study week). The timing of outbreak onset (On), peak (Pn.i), and resolution (Rn) are marked with red, green and purple bars, respectively, for the n th wave and the i th peak within that wave. If no resolution period occurs, the outbreak remains ongoing until reaching resolution or the next wave's onset (shown in Wave 3). Data used to develop this visualization are reported in Supplementary Excel Table S2. Critical periods included the acceleration period from onset to peak, the deceleration period from peak to resolution (or next outbreak onset if no resolution occurred), and the steady-state period from resolution to onset. Each period had a distinct behavior: acceleration periods indicated times when weekly rates steadily increased, deceleration periods indicated when rates steadily decreased, and steady-state periods indicated when rates plateaued. We also calculated the pace of increase (proxy to basic reproduction number) and pace of decline as the linear approximation from onset timing to peak timing and peak timing to resolution timing or next wave onset timing, respectively. We reported critical period duration in weeks and pace of increase or decline in cph per week (cph/week). By comparing the timing of critical points, duration of critical periods, and magnitude of smoothed peak rate times in different locations, we defined regional hotspots and explored the potential of traveling waves of infection outwards from regional hotspots.
To produce smoothed outbreak signatures and detect change points in the time series (including the time periods when the absolute change increased, decreased, or Middle panel: a heatmap that compactly reports weekly rates (outbreak magnitude) where light purple represents low values and dark purple represents high values. Bottom panel: a heatmap of the change in rates (trivial derivative) where acceleration periods are shown by the orange gradient, deceleration periods are shown by the blue gradient, and steady-state periods or critical change points near peaks are shown in grey. All panels share the common horizontal axis of time (study week). The timing of outbreak onset (O n ), peak (P n.i ), and resolution (R n ) are marked with red, green and purple bars, respectively, for the nth wave and the ith peak within that wave. If no resolution period occurs, the outbreak remains ongoing until reaching resolution or the next wave's onset (shown in Wave 3). Data used to develop this visualization are reported in Supplementary Excel Table S2. Critical periods included the acceleration period from onset to peak, the deceleration period from peak to resolution (or next outbreak onset if no resolution occurred), and the steady-state period from resolution to onset. Each period had a distinct behavior: acceleration periods indicated times when weekly rates steadily increased, deceleration periods indicated when rates steadily decreased, and steady-state periods indicated when rates plateaued. We also calculated the pace of increase (proxy to basic reproduction number) and pace of decline as the linear approximation from onset timing to peak timing and peak timing to resolution timing or next wave onset timing, respectively. We reported critical period duration in weeks and pace of increase or decline in cph per week (cph/week). By comparing the timing of critical points, duration of critical periods, and magnitude of smoothed peak rate times in different locations, we defined regional hotspots and explored the potential of traveling waves of infection outwards from regional hotspots.
To produce smoothed outbreak signatures and detect change points in the time series (including the time periods when the absolute change increased, decreased, or approached near-zero values), we applied a modified version of the Kolmogorov-Zurbenko (KZ) adaptive filter [37][38][39][40]. The KZ filter is a non-parametric smoothing technique used to detect distinct changes in a highly noisy time series [37][38][39][40]. This technique used a realvalues time series X(t), t = 0, ±1, ±2, ±3 . . . to generate a set of moving averages (MA) for windows of various sizes (Equation (3)): where q was the half length of a MA window, so Z(t, q) were the smoothed values of X(t) with a different degree of reduced noise. We created a set of time series of natural logarithm (ln)-transformed corrected rates (Equation (4)), estimated as: where C t,s were weekly counts for t-week and s-location corrected by a = 1 case and adjusted for estimated weekly population (P t,s ) estimates, per M = 100,000 persons. Natural logarithm transformations addressed the right skewed distribution of weekly rates due to short bursts of high values during epidemic peaks and prolonged periods of low reported rates in the beginning of the epidemic. We then created a set of time series of the trivial derivatives, which were estimated as the absolute change in weekly ln-transformed rates: s . We applied the KZ filter iteratively and generated a set of MA for both sets (Equations (5) and (6)), as: for q = 1-25, to cover various window sizes ω = 2q + 1 , from 3 to 51 weeks to depict the outbreak signature and to reflect the temporal changes from about one month to about one year. We then examined each set of select windows that sufficiently reduced noise yet maintained inflection points within the trivial derivative, as a verification step. By comparing the fit with sum of squares criteria, we selected candidates averaging across four windows: q = 1-15, q = 1-11, q = 1-7, and q = 3-5. We selected the best performing smoother (q = 3-5) to assess the similarities across governorate-specific outbreak signatures and to estimate the absolute change in weekly smoothed rates, as: Z t,s,3-5 = (Z * t+1,s,3 + Z * t+1,s,5 )/2 − (Z * t,s,3 + Z * t,s,5 )/2. The fit of this smoother for each governorate and nationally is shown in Figures S1-S6 (data in Supplementary Excel Table S3).
We estimated local maximums and minimums in Y * t,s , Z * t,s,3-5 , Z * t,s,0-7 , Z * t,s,0-11 , Z * t,s,0-15 , to identify time points and periods of high and low disease rates, respectively. To confirm maximums and minimums detected in the smoothed values of rates, we identified weeks when Y * t,s , Z * t,s,3-5 , Z * t,s,0-7 , Z * t,s,0-11 , and Z * t,s,0-15 approached zero. We defined outbreak onset as the beginning of a rapid rise of rates (Z * t,s,3-5 ) and inflection of the trivial derivative (Z t,s,3-5 ) from near-zero to high positive values. We defined outbreak peak as the time when rates were reaching a local maximum of Z * t,s,3-5 with the inflection point of Z t,s,3-5 at near-zero or negative values. We defined outbreak resolution time as when rates Z * t,s, [3][4][5] were low or near a local minimum and Z t,s,3-5 approached near-zero values. We expressed critical time points and periods in WHO-defined epidemiological weeks along with its plausible range.
We defined outbreak waves as either the duration from onset to resolution critical points or the duration between two onset critical points if no definite resolution occurred using the most refined window of 3-5 weeks. If no resolution occurred, we differentiated between waves by identifying weeks when rates Z * t,s,3-5 plateau for 4-6 weeks followed by 3-5 weeks of a steady increase in rates when the trivial derivative (Z t,s,3-5 ) increased from near-zero to high positive values. When possible, we identified primary and secondary outbreak peaks for specific waves based on our ability to detect local minimums and maximums. Secondary peaks referred to times when rates reached a local maximum of Z * t,s,3-5 after a burst of Z t,s,3-5 from high negative to high positive values. To confirm the potential for a secondary peak, we calculated the acceleration period and pace of increase from wave onset to each peak timing and calculated the deceleration period and pace of decline from each peak timing to the conclusion of the wave (resolution or following wave onset) that are characteristic for an ongoing outbreak.

Clustering of Outbreak Signatures
We applied pairwise Spearman cross-correlation estimates, covering six lead and six lag weeks (or lead-lag of +6 to −6 weeks) to examine similarities across governorate-specific outbreak signatures using the best performing smoother (Supplementary Table S1). We selected this lead-lag structure because governorate critical points varied by approximately six weeks. We assigned governorates with strong serial cross-correlations, close spatial proximity to each other, and similarities across outbreak signature features including the timing of critical points, duration of critical periods, and magnitude of peak rates into the same cluster (six clusters in all).
We performed all statistical analyses and visualizations using Excel (14.3.6), Stata (SE 15.1), and R (3.6.3) software. We created all maps using publicly available shapefiles reported by the United Nations (UN) Office for the Coordination of Humanitarian Affairs (OCHA) and published on the Humanitarian Data Exchange (HumData) [41]. The data and R codes used to create maps and illustrations in this study are available in Supplementary Materials.

Results
We presented results starting with quantifying outbreak features for Yemen nationally and by governorate using smoothed rates and trivial derivatives. We then examined similarities in outbreak features (critical points, critical periods, spatial proximity, and smoothed rate magnitudes) to define clusters. Finally, we analyzed similarities within and differences between clusters on a one-by-one basis in all subsequent subsections.

Four Waves of Cholera Outbreak Signatures
Over our 172 week study period, we identified~2.1 million confirmed cholera infections among the~28.6 million estimated inhabitants of Yemen. Based on the estimated onset, peak, and resolution critical points and acceleration, deceleration, and steady-state period durations for the Yemeni cholera outbreak from 4 September 2016 through 29 December 2019, we identified four distinct national outbreak waves ( Figure 2, Tables 1 and 2;  Supplementary Excel Tables S3 and S4). Wave 1 began in 12-18 September 2016 with a peak in 7-13 November 2016 and resolution in 2-8 January 2017. Peak rates were low during this wave (peak rate~0.06 cph) and the outbreak only affected 15 governorates. Wave 2, the largest outbreak wave nationwide, began approximately nine weeks after the conclusion of Wave 1 (6-12 March-2017). We identified two peaks in Wave 2: the first in 24-30 July 2017 (peak rate~151.3 cph) and the second in 18-24 September 2017 (peak rate~141.1), with a sharp decline in rates between peaks (~20% of peak rate values). The speed of increase for the primary and secondary peaks was 1.88-2.32 cph/week and 1.05-1.76 cph/week greater than the speed of decline, respectively. The prolonged decline from the second peak to the time of resolution (2-8 April 2018) lasted 28 weeks. For both Waves 1 and 2, the acceleration period from the onset to the peak was 1.5-2.2-fold shorter than the deceleration period from the peak to the onset of the next wave (Table 2). uary 2019 and Wave 4 from 28 January 2019 through 29 December 2019. In Wave 3, we identified two peaks: the first in 22-28 October 2018 (peak rate ~53.1 cph) and the second in 17-23 December 2018 (peak rate ~45.6 cph). These peaks were approximately 66% lower than the peak rate values for Wave 2, declined briefly (~10 cph) in rate values between peaks, and had a speed of increase/decline of ~2 cph/week. For Wave 3, the acceleration period from the onset to the peak was >1.5-fold longer than the deceleration period from the peak of Wave 3 to the onset of Wave 4. For Wave 4, we identified three peaks (22)(23)(24)(25)(26)(27)(28) April 2019, 15-21 July 2019, 16-22 September 2019). These peaks declined in magnitude throughout the outbreak wave (~99.7 cph, ~73.4 cph, ~66.1 cph) and occurred approximately 12 weeks apart. We did not identify a resolution critical point for Wave 4 by the end of our study.  Smoothed rates (laboratory-confirmed cholera cases per 100,000 persons) were estimated using the best performing smoother. Vertical bars represent the onset timing (red) for each wave (O1, O2, O3, O4), the peak timing (green) including one peak for Wave 1 (P1), two peaks for Wave 2 (P2.1, P2.2) and Wave 3 (P3.1, P3.2), and three peaks for Wave 4 (P4.1, P4.2, P4.3), and the resolution timing (purple) for Waves 1 and 2 (R1, R2). The heatmap of rates uses a purple gradient scale (lighter hues 0.0 cph and darker hues~400 cph) while the heatmap of derivative values uses a diverging scale ranging from −100 (dark blue) to +100 (dark orange) with values near zero in grey. All plots share a common horizontal axis of time reported in weeks since Week 36 of 2016 (defined as study Week 0). Governorates are grouped by assigned outbreak clusters: the core outbreak (Sana'a, Sana'a City, and Al-Hudaydah), immediate neighboring I (Amran, Al-Mahwit, Dhamar, and Al-Bayda), immediate neighboring II (Raymah, Ibb, and Taizz), northern (Hajjah, Al-Jawf, and Sa'ada), southern (Al-Dhale'e, Abyan, Aden, and Lahj), and eastern (Marib, Al-Maharah, and Shabwah) clusters. Data used to develop this visualization are reported in Tables 1 and 2 and Supplementary Excel Tables S3 and S4. All critical points were estimated using smoothed rates (laboratory-confirmed cholera cases per 100,000 persons) and trivial derivative values (∆ rates per epidemiological week) produced with Kolmogorov-Zurbenko (KZ) adaptive filters. We selected the best performing smoother (Z * t,s,3-5 ) as the average of a 3 week (Z * t+1,s,3 ) and 5 week (Z * t+1,s,5 ) window size and trivial derivative as the absolute change in weekly smoothed rates calculated as: Z t,s, [3][4][5]    Deceleration periods are reported for both the duration from a wave's peak to resolution and from a wave's peak to the following wave's onset. Acceleration periods indicated times when weekly rates were steadily increasing (onset to peak), deceleration periods indicated times when weekly rates were steadily decreasing (peak to resolution or peak to onset), and steady-state periods indicated times when the weekly rates plateaued (resolution to onset). The pace of increase or decline is measured as a linear slope (rates per week) and estimated simply as the smoothed rate at one critical point minus the smoothed rate at the following critical point and divided by the duration between critical points. Duration between critical points ranged 2-56 weeks. Critical periods are denoted from their starting and ending onset (O n.i ), peak (P n.i ), and resolution (R n ) critical points for the nth outbreak wave and ith occurrence per wave. No resolution critical point for Wave 4 was found by the conclusion of the study period. As a result, we found no deceleration or steady-state periods for Wave 4. Critical period duration is reported in weeks since Week 36 of 2016 (defined as study Week 0). The pace of increase and decline are reported in laboratory-confirmed cholera cases per 100,000 persons per week. Critical points used to define critical periods were estimated using smoothed rates (laboratory-confirmed cholera cases/100,000 persons) and trivial derivative values (∆ rates/week) produced with Kolmogorov-Zurbenko (KZ) adaptive filters (see Table 1). Governorates are grouped by assigned outbreak clusters indicated with double solid lines and include: the core outbreak (Sana'a, Sana'a City, and Al-Hudaydah), immediate neighboring I (Amran, Al-Mahwit, Dhamar, and Al-Bayda), immediate neighboring II (Raymah, Ibb, and Taizz), northern (Hajjah, Al-Jawf, and Sa'ada), southern (Al-Dhale'e, Abyan, Aden, and Lahj), and eastern (Marib, Al-Maharah, and Shabwah) clusters. Data summarized in this table are, in part, graphically presented in Figure 2 and Supplementary Figures S1-S6.
We also identified two additional waves: Wave 3 from 28 May 2018 through 27 January 2019 and Wave 4 from 28 January 2019 through 29 December 2019. In Wave 3, we identified two peaks: the first in 22-28 October 2018 (peak rate~53.1 cph) and the second in 17-23 December 2018 (peak rate~45.6 cph). These peaks were approximately 66% lower than the peak rate values for Wave 2, declined briefly (~10 cph) in rate values between peaks, and had a speed of increase/decline of~2 cph/week. For Wave 3, the acceleration period from the onset to the peak was >1.5-fold longer than the deceleration period from the peak of  We also found variation in the number and timing of peaks for Wave 3 (12 governorates with one peak, 6 governorates with two peaks) and Wave 4 (1 governorate with one peak, 8 governorates with two peaks, and 11 governorates with three peaks).

Six Clusters of Outbreak Signatures
Based on the combination of correlation estimates, timing of critical points, and geographic proximity, we assigned all governorates into one of six clusters including a core outbreak cluster, two immediately neighboring clusters, and three remote clusters in northern, southern and eastern regions of Yemen (Figure 3; Supplementary Excel Table S3). Smoothed rates (laboratory-confirmed cholera cases per 100,000 persons) were estimated using the best performing smoother. Vertical red, green, and purple lines depict the timing of onset, peak, and resolution critical point (s), respectively, for each of the four outbreak waves. Bottom panels: governorate-level signatures are shown in two columns of stacked multi-panel time series plots that share a common vertical axis between columns (log10-transformed rates (cph)) and a common horizontal axis within columns (weeks since Week 36 of 2016, defined as study Week 0). Vertical green and purple bars depict the timing of the Wave 2 peak, Wave 2 resolution, and Wave 3 peak, which were important signature features for defining outbreak clusters. Spearman cross correlation estimates (ρ) are at lag 0 between adjacently reported governorates within the same cluster based on 172 weeks (all correlation coefficients are significant at α < 0.001). Data used to develop this visualization are reported in Table 1, Supplementary Table S1, and Supplementary Excel Table S3.

Core Cluster
We identified Sana'a, Sana'a City, and Al-Hudaydah governorates as the core outbreak cluster (Figure 4; Supplementary Excel Table S5). We used the term 'core' as these governorates were located centrally in Yemen, contained or surrounded the nation's capital, housed a large percentage of the population, and had the earliest timing of onset and peak critical points in Wave 1. These governorates had closely aligned primary and sec-  Smoothed rates (laboratory-confirmed cholera cases per 100,000 persons) were estimated using the best performing smoother. Vertical red, green, and purple lines depict the timing of onset, peak, and resolution critical point (s), respectively, for each of the four outbreak waves. Bottom panels: governorate-level signatures are shown in two columns of stacked multi-panel time series plots that share a common vertical axis between columns (log10-transformed rates (cph)) and a common horizontal axis within columns (weeks since Week 36 of 2016, defined as study Week 0). Vertical green and purple bars depict the timing of the Wave 2 peak, Wave 2 resolution, and Wave 3 peak, which were important signature features for defining outbreak clusters. Spearman cross correlation estimates (ρ) are at lag 0 between adjacently reported governorates within the same cluster based on 172 weeks (all correlation coefficients are significant at α < 0.001). Data used to develop this visualization are reported in Table 1, Supplementary Table S1, and Supplementary Excel Table S3.

Core Cluster
We identified Sana'a, Sana'a City, and Al-Hudaydah governorates as the core outbreak cluster (Figure 4; Supplementary Excel Table S5). We used the term 'core' as these governorates were located centrally in Yemen, contained or surrounded the nation's capital, housed a large percentage of the population, and had the earliest timing of onset and peak critical points in Wave 1. These governorates had closely aligned primary and secondary peaks in Wave 2 from 26 June 2017 to 30 July 2017 (131.4-172.8 cph) and from 18 September 2017 to 1 October 2017 (82.1-215.8 cph). Pace of increase for primary peaks were 2.43-2.75-fold greater than pace of decline for all governorates. After a long deceleration period of~30 weeks, cholera rates never returned to pre-Wave 2 baseline levels for any of these three governorates. Rates decreased by~40-66% for these governorates during the~6 weeks between peaks. While Sana'a had a smaller third peak (66.8 cph) in 1-7 January 2018, Sana'a City and Al-Hudaydah had a brief (9-14 weeks) resolution period from 19 March 2018 to 1 July 2018.   Table 1 and Supplementary Excel  Table S5.  Al-Hudaydah:~81 cph). In Wave 4, these governorates had the same first peak timing in 15-28 April 2019. However, the second peak in Sana'a and Sana'a City occurred in 22-28 July 2019, approximately two months prior to the third peak of Sana'a City and second peak of Al-Hudaydah in 16-22 September 2019. Differences in the alignment of peak timing supported Sana'a and Sana'a City's strong association (ρ = 0.82-0.89; p < 0.001) with Al-Hudaydah for lags from 1 to 5 weeks.

Immediate Neighboring Cluster I
We identified four governorates, Amran, Al-Mahwit, Dhamar, and Al-Bayda, that formed a cluster immediately neighboring the core cluster ( Figure 5; Supplementary Excel Table S5). This neighboring cluster had similar critical points but higher peak rates in Wave 2 compared to the core cluster. In Wave 3, all three governorates had two peaks: from 24 September 2018 to 21 October 2018 and from 19 November 2018 to 30 December 2018. For each governorate, primary and secondary peak rate values were nearly identical (Sana'a: ~112 cph; Sana'a City: ~51 cph; Al-Hudaydah: ~81 cph). In Wave 4, these governorates had the same first peak timing in 15-28 April 2019. However, the second peak in Sana'a and Sana'a City occurred in 22-28 July 2019, approximately two months prior to the third peak of Sana'a City and second peak of Al-Hudaydah in 16-22 September 2019. Differences in the alignment of peak timing supported Sana'a and Sana'a City's strong association (ρ = 0.82-0.89; p < 0.001) with Al-Hudaydah for lags from 1 to 5 weeks.

Immediate Neighboring Cluster I
We identified four governorates, Amran, Al-Mahwit, Dhamar, and Al-Bayda, that formed a cluster immediately neighboring the core cluster ( Figure 5; Supplementary Excel Table S5). This neighboring cluster had similar critical points but higher peak rates in Wave 2 compared to the core cluster.   Table 1 and Supplementary Excel Table S5.
Though the early primary peak (24 July 2017 to 27 August 2017) was similar to the core cluster, the neighboring cluster had a later secondary peak from 5 March 2018 to 1 April 2018 that fell~10 weeks after the core cluster. Additionally, Al-Mahwit, Al-Bayda, and Amran had peak rates of~400 cph, roughly twice the magnitude of peak rates in the core cluster. The pace of increase for primary peaks ranged from 10.8 to 20.9 cph/week, which was on average twice as large compared to the core cluster (6.6-10.8 cph/week). Unlike the core cluster, none of the governorates in this neighboring cluster had a Wave 2 resolution, though rate values declined to~20 cph (roughly twice the magnitude of the core cluster) by the onset timing of Wave 3.
In Wave 3, Amran, Al-Mahwit, and Dhamar had only one peak, which occurred from 24 September 2018 to 21 October 2018. Al-Bayda had a secondary peak in 19-25 November 2018, which was identical in timing to governorates in the core cluster. As in Wave 2, Wave 3 peak rate values for the neighboring cluster were, on average, higher than the peak rates for core cluster governorates (63.1-188.4 cph).
In Wave 4, all governorates peaked twice (except Al-Bayda, which peaked three times). All four governorates shared a secondary peak in 16-29 September 2019, which aligns with the later peak timing for Sana'a City and Al-Hudaydah.

Immediate Neighboring Cluster II
Raymah, Ibb, and Taizz formed a second cluster that immediately neighbored the core cluster. This second neighboring cluster had two well-aligned peaks in Wave 2, one wellaligned peak in Wave 3, three peaks in Wave 4, and low peak rate magnitudes compared to previously described clusters (Figure 6; Supplementary Excel Table S5).
In Wave 2, the primary peak occurred from 26 June 2017 to 23 July 2017 (like the core cluster) and the secondary peak occurred from 15 January 2018 to 11 March 2018 (like the first neighboring cluster). Wave 2 peak rates ranged from 9.25 to 125.6 cph, which were approximately half the magnitude of the first neighboring cluster (216.3-397.7 cph). All governorates in this second neighboring cluster reached resolution timing for Wave 2 in 2-15 April 2018, when rates nearly returned to pre-Wave 2 levels (~2.45-3.56 cph).
In Wave 3, the governorates in this neighboring cluster had only one outbreak peak, which occurred from 24 September 2018 to 25 November 2018 (same as core cluster). However, peak rates had low magnitudes, ranging from 22.3 to 47.2 cph, which were less than half the peak rate magnitudes for governorates in the core and first neighboring clusters. In Wave 4, this neighboring cluster had similar peak timing as the core cluster for all three identified peaks.

Remote Clusters: Northern, Southern, and Eastern
The ten remaining governorates had outbreak signatures that differed from the governorates in the core and neighboring clusters, forming three geographically distinct clusters in the northern (Hajjah, Al-Jawf, and Sa'ada), southern (Al-Dhale'e, Abyan, Aden, and Lahj), and eastern (Marib, Al-Maharah, and Shabwah) regions of Yemen. Governorates within these remote clusters had one peak in Waves 2 and 3, a short Wave 2 deceleration period duration, and an early arrival of Wave 2 resolution critical points compared to the other three clusters.
In the northern cluster (Figure 7; Supplementary Excel Table S5), governorates had Wave 2 peak timing from 4 September 2017 to 8 October 2017. Rate values reached magnitudes like those in the core outbreak cluster (73.0-208.1 cph). However, the pace of increase to and pace of decline from Wave 2 peak timing were approximately the same and oc-curred~6-10 weeks after the primary peak timing of governorates in the first neighboring cluster. Cholera rates declined to pre-Wave 2 magnitudes (<1 cph) 16-30 weeks after the primary peak. Governorates within this northern cluster had Wave 3 peak rate magnitudes (24.0-49.0 cph) like those in the second neighboring cluster governorates.  Table 1 and Supplementary Excel Table S5.  Table 1 and Supplementary Excel Table S5.
Wave 2 peak timing from 4 September 2017 to 8 October 2017. Rate values reached magnitudes like those in the core outbreak cluster (73.0-208.1 cph). However, the pace of increase to and pace of decline from Wave 2 peak timing were approximately the same and occurred ~6-10 weeks after the primary peak timing of governorates in the first neighboring cluster. Cholera rates declined to pre-Wave 2 magnitudes (<1 cph) 16-30 weeks after the primary peak. Governorates within this northern cluster had Wave 3 peak rate magnitudes (24.0-49.0 cph) like those in the second neighboring cluster governorates.   Table 1 and Supplementary Excel Table S5. In the southern cluster ( In the southern cluster (Figure 8; Supplementary Excel Table S5), Al-Dhale'e, Abyan, and Aden had only one Wave 2 peak that occurred in 17-23 July 2017. Lahj had two Wave 2 peaks: the primary occurred in 14-20 August 2017 and the secondary occurred in 25-31 December 2017. Irrespective of peak timing, all governorates had a prolonged deceleration period of ~30 weeks and reached Wave 2 resolution from 22 January 2018 to 18 March 2018. By the conclusion of Wave 2, governorates within this cluster had returned to pre-Wave 2 rate magnitudes (<1 cph).   Table 1 and Supplementary Excel Table S5.
Governorates in the eastern cluster (Figure 9; Supplementary Excel Table S5) had no reported cases of cholera during Wave 1 and no well-defined critical points in Wave 3. In fact, cholera rates did not exceed 5 cph between the resolution of Wave 2 (as early as [18][19][20][21][22][23][24] December 2017) and the onset of Wave 4 (as late as 18-24 March 2019) for any eastern cluster governorate. In Wave 2, these governorates had a wide range of peak rate magnitudes, from Shabwah at~22.1 cph to Marib at~110.2 cph. Shabwah and Al-Maharah had the earliest resolution timing, which occurred from 18 December 2017 to 4 February 2018. As with governorates in the southern cluster, all governorates in the eastern cluster had Wave 4 peak rate magnitudes of <35 cph.  Table 1 and Supplementary Excel Table S5.  Table 1 and Supplementary Excel Table S5.

Discussion
Our study demonstrated the utility of publicly reported surveillance data to characterize, classify, and compare infectious disease outbreak signatures. We proposed an application of non-parametric KZ adaptive filters to define outbreak signatures and identified regional hotspots of cholera outbreaks. We used this data-driven approach to define outbreak critical points (onset, peak timing, and resolution) and critical periods (acceleration, deceleration, and steady state) to examine spatiotemporal patterns in cholera outbreak signatures. This methodology introduced a protocol for routine vulnerability mapping of outbreak hotspots to improve resource management and mobilization during humanitarian aid responses.
Application of KZ filter methodology permitted closer inspection of subtle differences between governorates' outbreak signatures considering several investigated features. Using the proposed methodology, we defined a core outbreak hotspot of Sana'a, Sana'a City, and Al-Hudaydah, which was the expected origin of the national outbreak [42][43][44][45]. These governorates shared nearly identical features as the national outbreak signature in all four waves, including onset and peak timing critical points, peak rate magnitudes, and durations of acceleration and deceleration periods. While key informant interviews and trend analyses elsewhere suggested that Wave 1 peak timing occurred in 5-31 December 2016, we identified earlier peak timing for the core cluster and nationally ranging from 10 October 2016 to 13 November 2016 [26,27,46]. Additionally, various sources including the WHO declared that Wave 1 concluded for all governorates and nationally in 13-27 April 2017 [26,27,46,47]. However, we identified a resolution period for all affected governorates ranging from 2 January 2017 to 5 March 2017.
Our results for Wave 2 outbreak signatures closely aligned with those of Dureab et al. and Camacho et al., whose descriptive trend analyses were based on the Yemeni National Electronic Disease Early Warning System (eDEWS) [26,27]. While Camacho et al.  [26,27]. We identified two outbreak peaks with similar timing (24-30 July 2017 and 25 September 2017 to 01 October 2017), shared by most governorates. However, we also identified a third peak found in the two neighboring clusters (Amran, Al-Mahwit, Dhamar, Al-Bayda, Raymah, Ibb, and Taizz) surrounding the core outbreak cluster (Sana'a, Sana'a City, and Al-Hudaydah) from 5 March 2018 to 1 April 2018. Similarities in our findings suggested consistency in the reporting of district-level and governorate-level data.
These findings potentially suggested traveling waves from the core outbreak cluster, whose earlier second peak timing (18 September 2017 to 1 October 2017) preceded the later peak timing (1 January 2018 to 20 May 2018) in neighboring governorates. The core cluster reflected a persistent cluster of infections due to consistently high peak rates, early peak timing, and large population density in these governorates. The later peak timing of governorates within close spatial proximity suggested a high degree of percolation to other clusters, though these variations may have been dictated by differences in the spatial distribution of conflict, environmental, and socioeconomic factors within each governorate. Further research is needed to describe patterns of disease spread and differences in drivers of infection in each governorate.
Our findings also suggest a possible second traveling wave from the first neighboring cluster (Amran, Al-Mahwit, Dhamar, and Al-Bayda) to the northern, southern, and eastern remote clusters. The four governorates in this neighboring cluster bordered governorates in each of the three remote clusters. As described earlier, this neighboring cluster is characterized by extremely high peak rate magnitudes (216.3-397.7 cph) for the first peak timing (24 July 2017 to 27 August 2017) in Wave 2 with equally high paces of increase to the peak (10.8-20.9 cph/week). In contrast, most remote clusters had a single peak whose timing was slightly later than the neighboring cluster (24 July 2017 to 17 September 2017) and whose magnitude was lower (22.14-164.88 cph) in most governorates. This suggested that larger outbreaks in this neighboring cluster could have resulted in percolation of the outbreak to the remote clusters, which experienced a short, intense Wave 2 outbreak before quickly returning to pre-outbreak rate values (<4 cph) from 18 December 2017 to 1 April 2018.
While WHO epidemiological surveillance bulletins suggested that Wave 2 persisted from April 2017 through January 2020, we identified two national outbreak waves with onsets in May 2018 and January 2019 [47]. These waves suggest that the epidemic remains ongoing with numerous onsets of outbreak events rather than one prolonged outbreak period. For Waves 1 and 2, the acceleration period from the onset to the peak was 1.5-2.2-fold shorter than the deceleration period from the peak timing to the onset of the next wave. Furthermore, for Wave 3, the acceleration period from the onset to peak was >1.5-fold longer than the deceleration period from the peak of Wave 3 to the onset of Wave 4. As with Waves 1 and 2, we found great heterogeneity in the onset, peak timing, and duration of outbreaks across governorates for Wave 3 and Wave 4. Differences in patterns of outbreak signatures across waves may suggest that drivers of disease transmission varied over time and location. These differences must be inspected individually for each wave and may have required different types of humanitarian relief according to the conflict-, environmental-, and socioeconomic-related drivers of infection.
Discrepancies in outbreak signatures can significantly influence reported associations between cholera rates and human-made or environmental drivers of transmission. For example, Wave 3 began nationally from 28 May 2018 to 3 June 2018, though onset timing ranged from 2 April 2018 to 8 July 2018. The Wave 3 onset timing for southern cluster governorates closely aligned with the arrival of the Mekunu cyclone, which made landfall in Oman and Yemen on 25 May 2018 and resulted in severe flooding [48][49][50]. In contrast, the Battle of Al-Hudaydah from 13 to 22 June 2018 left governorates in the core, neighboring, and remote northern clusters with severe health infrastructure destruction [51]. These natural disaster and conflict events, whose repercussions persisted through December 2019, may explain the shorter acceleration periods and more frequent outbreak peaks in Waves 3 and 4. Future research using our defined outbreak signatures is needed to examine how environmental and human-made risk factors varied both across and within outbreak clusters on a wave-by-wave basis.
Our study could shed light on previous research examining associations between cholera incidence and environmental risk factors for Wave 2 of the Yemeni cholera outbreak. Camacho et al. used district-level surveillance records and high-resolution rainfall data to show that within 10 days of a 10 mm rainfall event the likelihood of cholera infection increased by 21% compared to a week of no rainfall [26]. However, results assumed that each district followed identical outbreak onset and peak timing according to nationallevel, WHO-defined outbreak wave critical points [26]. Other studies aimed at modeling cholera outbreak signatures similarly did so with pre-defined national-level characteristics, preventing exploration of regional hotspots [27,46,52]. Our results provided more refined estimates of critical points and helped account for the variability of outbreak signatures or differences in outbreak duration across governorates, which faced differing intensities of wartime conflicts, food insecurity, and malnutrition. If more spatially granular time series data were reported at the district level, we could further refine how governorates were assigned to outbreak clusters.
Our study had several limitations. Firstly, the reliability of case data may be impacted by the scattered dispersal of public technical reports and their inconsistent, infrequent, and sometimes non-overlapping time periods. Furthermore, data collection and reporting were expected to be limited by health care services/capacity, utilization of health services, and reporting ability in combination with the ruralness and governing stability of each governorate. The WHO provided no metadata explaining limitations in data quality and coverage after data were aggregated from district-level eDEWS reports to governorate-level WHO bulletins. Dureab et al., who used eDEWS data reports, found only~1 week lag between the time of infection and reporting of time series records suggesting consistent reporting despite ongoing conflict events [27]. However, key informant interviews in Yemen noted significant destruction of emergency health facilities in Sana'a City and low capacity for routine and reliable laboratory culture testing in Sana'a and Al-Hudaydah [3,5,46]. This suggests possible under-reporting in these governorates, especially during the second wave of the epidemic between March 2017 and April 2018.
We compiled data presented in various formats to the best of our knowledge and abilities and are unaware of alternative methods that could have better validated the compiled data, given the inconsistent reporting and availability of case reports. Similarly, data were only available at the governorate level, which prevented more refined cluster assignments that would account for district-level differences in topography, climate, health infrastructure, etc. While district-level data would improve the granularity of classified outbreak clusters, our use of governorate-level data reflected both the usability of publicly reported data and the reality of the data landscape in Yemen. That said, we found similar critical point timing and magnitudes in our national outbreak signature to the study performed by Dureab et al. (which used eDEWS data), which suggested that our estimation techniques retained the features of reported data [27].
Secondly, population estimates used to calculate weekly rates only accounted for possible changes in population growth rate and conflict fatalities. This growth rate was retrieved from the 2004 Yemeni Census as no additional population estimates or growth rates were provided in the 2014 Census [35,53]. ACLED fatality data used to adjust for violent and non-violent political events were retrieved from a variety of media sources and may not be reliable [54][55][56][57]. Additionally, we were unable to adjust for possible population migration from Yemen or internal displacement within Yemen, as these data were not reported consistently at the monthly or weekly levels [36]. While a substantial proportion (~15%) of Yemen's population has been internally displaced throughout the Yemeni Civil War, neither ACLED nor any other dataset we explored provided data on migration or internal displacement with sufficient temporal resolution, spatial granularity, or consistent reporting to use in adjusting population estimates. Furthermore, data lacked sufficient spatial resolution to consider examining rates by population density, as most governorates included both urban cities and rural mountainous areas. Because we used a prorated weekly population estimate to account for possible population growth over time in the absence of more refined population data, we likely underestimated the rate of infection.
Finally, data from this study were reported through 29 December 2019. We were unable to extend our time series through 2020 or 2021 due to a lack of epidemiological bulletins for this year [29,47]. Understandably, resources could have been diverted to mitigate and manage the ongoing SARS-CoV-2 pandemic. While weekly updates were provided on the current cholera crisis in Yemen, these reports provided insufficient temporal and spatial granularity for recreating weekly time series of cases. However, various reports suggested, as found in our study, that the cholera epidemic within Yemen remains ongoing with rates like those reported in 2019 [58][59][60][61][62]. Thus, coordinated efforts for mitigating and preventing the spread of cholera are still needed [62][63][64][65][66].
While we detected potential traveling waves for governate-specific outbreaks, there could be many reasons for observed spatiotemporal patterns. The differences in outbreak signatures could be governed by variations in environmental factors, such as rainfall, flooding, and water contamination, conflicts and social unrest, depleted health infrastructure, food insecurity, and many other factors that deserve further investigations. Any such investigation should start with better understanding the spatiotemporal behaviors of infection across clusters and synchronicity between exposure events and outbreaks, that could be explored with dynamic mapping. Dynamic or animated maps are an emerging data visualization technique that compresses large amounts of spatially and temporally aligned records into visuals. We have demonstrated the utility of this visualization technique to identify persistent and percolating clusters of Salmonella outbreaks and traveling waves of influenza in the United States [18,67]. Our curated dataset provides the foundation for aligning and integrating other spatiotemporal data on conflict-, environmental-, and socioeconomic-related factors to conduct these future analyses. This expanded dataset can then be used for evaluating the drivers and consequences of disease outbreaks during the Yemeni Civil War.
To increase the usability and transparency of valuable epidemiological records collected by many national and international organizations, public access and dissemination of information are critical. Publicly available data and information could assist in preventing or mitigating humanitarian crises such as this one in Yemen. We encourage the holders of epidemiological records to share data and metadata with the broader scientific community at the highest temporal and spatial granularity possible to ease the extraction, understanding, and utilization of publicly disseminated data. We also encourage future applications of the novel methods for early outbreak detection, regional hotspot identification, and resource distribution coordination in accordance with the global mandate of the WHO's GTFCC [8][9][10][11]. Our findings, in conjunction with future research, could support the GTFCC roadmap and provide solid methodology to define outbreak signatures, classify regional hotspots, evaluate risk factors contributing to outbreak transmission, and establish an early outbreak warning system [10,11].

Conclusions
Our study provided a data-driven approach to describe and compare outbreak signatures of governorate-level cholera outbreaks in Yemen. These techniques utilized time series surveillance data and were not unique to the temporal resolution, spatial granularity, or infectious disease examined. Though our analysis evaluated data retrospectively, the proposed methodology can be adapted and applied for real-time forecasting to determine outbreak signatures of ongoing epidemics. Furthermore, our approach demonstrated how data published publicly by humanitarian organizations can be retrieved, harmonized, and utilized for modeling outbreak signatures. The detected heterogeneity in outbreak signatures associated with spatial proximity suggested possible traveling waves of infection outwards from the core cluster. Further research is needed to investigate if these differences are due to impacts of conflict, environmental, social, and economic factors across governorates. Similarities in outbreak signatures also suggested potentially shared environmental and human-made drivers of infection. By identifying and grouping governorates with similar critical points, periods, and rate magnitudes, we provided the context for identifying and comparing human-made and environmental drivers of disease transmission. These methods can help to parse out associations between infectious disease outbreaks and wartime conflict, rainfall, food insecurity, and consumer price risk factors to assist in coordinating humanitarian relief efforts.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijerph19010378/s1, Supplementary Excel File: This file provides Excel Tables S1-S5, which publicly report datasets used for statistical analyses and visualizations within this study; Figures S1-S6: Stacked multi-panel time series plots of the national Yemeni cholera outbreak signature and the governorate-level signatures of each outbreak cluster; Table S1: